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Abstract 

In this paper, we propose a general approach called Generalized Multiscale Finite Element Method 
(GMsFEM) for performing multiscale simulations for problems without scale separation over a complex 
input space. As in multiscale finite element methods (MsFEMs) , the main idea of the proposed approach 
is to construct a small dimensional local solution space that can be used to generate efficient and accurate 
approximation to the multiscale solution with a potentially high dimensional input parameter space. In 
the proposed approach, we present a general procedure to construct the offline space that is used for a 
systematic enrichment of the coarse solution space in the online stage. The enrichment in the online stage 
is performed based on a spectral decomposition of the snapshot space consisting of all plausible functions 
in the space. In the online stage, for any input parameter, a multiscale space is constructed to solve 
the global problem on a coarse grid. The online space is constructed via a spectral decomposition of the 
offline space and by choosing the eigenvectors corresponding to the largest eigenvalues. The computational 
saving is due to the fact that the construction of the online multiscale space for any input parameter is fast 
and this space can be re-used for solving the forward problem with any forcing and boundary condition. 
Compared with the other approaches where global snapshots are used, the local approach that we present 
in this paper allows us to eliminate unnecessary degrees of freedom on a coarse-grid level. We present 
various examples in the paper and some numerical results to demonstrate the effectiveness of our method. 

1 Introduction 

1.1 Multiscale problems, the input-output relation, and the need for model 
reduction 

Many problems arising from various physical and engineering applications are multiscale in nature. Because 
of the presence of small scales and uncertainties in these problems, the direct simulations are prohibitively 
expensive. Moreover, these problems are typically solved for many source terms with input parameter 
coming from a high dimensional parameter space. For example, the flow in heterogeneous porous media 
described by Darcy's equation is typically solved for multiple source terms. Moreover, the permeability 
usually has uncertainties which are parametrized in some sophisticated manner. In this case, one needs 
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Figure 1: Flow chart 



to solve many forward problems with different source terms and a wide range of permeabilities to make 
accurate predictions. These problems can be cast using an input-output relation (see Figure Q} which is 
typically done in reduced-order modeling. For the example of flow problems, the input space consists of 
source terms and the permeability that takes a value from a large parameter space. The output space 
depends on the quantities of interest and may consist of coarse-grid solutions or some other integrated 
quantities with respect to the solution. In many applications, the output space is typically smaller than 
the input space. The design of a general multiscale finite element framework that takes advantage of the 
effective low dimensional solution space for multiscale problems with high dimensional input space is the 
main objective of this paper. 

Due to a large number of forward simulations, the computational effort can be tremendous to learn 
and process the output space given the high dimensionality of the input parameter space. In many of 
these problems, the solution space can be approximated by a low dimensional manifold via some model 
reduction tools. The main objective of reduced-order models is to represent the solution space with a 
small dimensional space. However, many existing reduced-order methods fail to give a small dimensional 
output solution space when the physical solution has multiscale structures. Another major limitation of 
the current reduced-order methods is that the reduced solution space needs to be regenerated for different 
forcing or boundary conditions. The general multiscale finite element framework proposed in this paper is 
designed to remove these two limitations by dividing the construction of reduced basis into the offline and 
online steps, and constructing our online multiscale bases from a reduced localized offline solution space. 

1.2 Local and global model reduction concepts 

Many local, global, and local-global model reduction techniques have been developed. The main idea of 
these methods is to find a small dimensional space that can represent the solution space given the input 
space. 

Global model reduction techniques (see e.g., [53l EH [HI IHJ [6j [55]) construct a space of global fields 
that can approximate the solution space. One can, for example, consider a space of exhaustive global 
snapshots obtained by solving the global problem for many input parameters. This space can be further 
reduced using a spectral decomposition. In practice, the resulting space is constructed by solving global 
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Figure 2: Flow chart 



problems for some selected input parameters, right hand sides, and boundary conditions. These methods 
have been used with some success in practice. However, when the right hand sides or boundary conditions 
are changed, the resulting reduced space must be recomputed. 

Local approaches (e.g., see [351 111 ID ED SI 13 ED H3 US EH ED] for upscaling and multiscale methods) at- 
tempt to approximate the solution in local (coarse-grid) regions for all input parameters without computing 
global snapshots of solutions. Local approaches first compute an offline space (possibly small dimensional) 
which is used to compute multiscale basis functions at the online stage. The local approximation space at 
the online stage is computed by finding a subspace of offline space for a given input parameter (see Figure 



Local approaches can be effective as they avoid the computation of global snapshots. Local approaches 
become more effective if the restriction of the solution space onto a local region has a small dimension. 
This is the case if the dimension of the space of solutions restricted to a coarse region is smaller than the 
dimension of the fine-grid space within this coarse region. For example, if the parameter is a coarse-grid 
scalar function, then at the coarse-grid level, this parameter is a scalar. While if we consider this problem 
from the point of view of a global model reduction, then the parameter belongs to a large dimensional 
space and this may not be amenable to computations. 

One of advantages of local approaches is that they eliminate the unnecessary degrees of freedom in the 
parameter space at the coarse-grid level. In global methods, one first needs to compute many expensive 
global snapshots and many snapshots may not contribute to the solution at the online stage. In local 
approaches, these values of the parameter are identified at the coarse level inexpensively. Moreover, local 
approaches can easily handle large-scale parameter space when the parameter is a coarse-grid function and 
local approximation spaces are usually independent of the source terms or boundary conditions. We will 
further elaborate these issues in the paper. 
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1.3 This paper 

In this paper, we introduce a general multiscale framework, which we call the Generalized Multiscale Finite 
Element Method (GMsFEM). This method incorporates complex input space information and the input- 
output relation. It systematically enriches the coarse space through our local construction. Our approach, 
as in many multiscale and model reduction techniques, divides the computation into two stages: offline; 
and online. In the offline stage, we construct a small dimensional space that can be efficiently used in the 
online stage to construct multiscale basis functions. These multiscale basis functions can be re-used for any 
input parameter to solve the problem on a coarse-grid. Thus, this provides a substantial computational 
saving at the online stage. Below, we present an outline of the algorithm and a chart that depicts our 
algorithm in Figure [5J 

1. Offline computation: 

— 1.0. Coarse grid generation; 

— 1.1. Construction of snapshot space that will be used to compute an offline space. 

— 1.2. Construction of a small dimensional offline space by performing dimension reduction in the 
space of global snapshots. 

2. Online computations: 

— 2.1. For each input parameter, compute multiscale basis functions; 

— 2.2. Solution of a coarse- grid problem for any force term and boundary condition; 

— 2.3. Iterative solvers, if needed. 

In the offline computation, we first set up a coarse grid where each coarse-grid block consists of a 
connected union of fine-grid blocks. The construction of snapshot space in Step 1.1 involves solving local 
problems for various choices of input parameters. This space is used to construct the offline space in Step 
1.2 via a spectral decomposition of the snapshot space. The snapshot space in a coarse region can be 
replaced by the fine-grid space associated with this coarse space; however, in many applications, one can 
judiciously choose the space of snapshots to avoid expensive offline space construction. The offline space 
in Step 1.2. is constructed by spectrally decomposing the space of snapshots. This spectral decomposition 
is typically based on the offline eigenvalue problem. The spectral decomposition enables us to select the 
high-energy elements from the offline space by choosing those eigenvectors corresponding to the largest 
eigenvalues. More precisely, we seek a subspace of the snapshot space such that it can approximate any 
element of the snapshot space in the appropriate sense defined via auxiliary bilinear forms. 

In the online step 2.1, for a given input parameter, we compute the required online coarse space. In 
general, we want this to be a small dimensional subspace of the offline space. This space is computed 
by performing a spectral decomposition in the offline space via an eigenvalue problem. Furthermore, the 
eigenvectors corresponding to the largest eigenvalues are identified and used to form the online coarse 
space. The online coarse space is used within the finite element framework to solve the original global 
problem. Here, we propose several options such as the Galerkin coupling of multiscale basis functions, 
the Petrov- Galerkin coupling of multiscale basis functions, etc. In some of these coupling approaches, the 
choice of the initial partition of unity (that can be computed in the offline or online stage) is important 
and it will be discussed in the paper. 

Our techniques differ from many previous approaches that are based on the homogenization theory. In 
the homogenization based methods, one usually constructs local approximation based on local solves and 
these approaches do not provide a systematic procedure to complement the local spaces. It is important 



4 



to note that one needs to systematically complement the local spaces in order to converge to the fine-grid 
solution. How to develop an online systematic enrichment procedure and how to construct the initial 
partition of unity functions play a crucial role in obtaining a low dimensional offline space. These issues 
are central points of our proposed method. 

We also discuss iterative solvers that use the coarse spaces and iterate on the residual to converge to 
the fine-scale solution. These iterative solvers serve as an online correction of the coarse-grid solution. We 
consider two-level domain decomposition preconditioners and some other methods where the importance 
of appropriately chosen multiscale coarse spaces has been demonstrated in the literature [531 ES] • We will 
discuss how the choice of coarse spaces yields optimal iterative solvers where the number of iterations is 
independent of the high contrast in the media properties. 

2 A Generalized Multiscale Finite Element Method. 

To describe the GMsFEM for linear problems, we consider 

= f, (1) 
subject to some boundary conditions, where /x is the parameter. For example, 

Ln(u) = — cHv(k(x; |ti)Vu). (2) 

Here, the operator L may depend on various spatial fields, e.g., heterogeneous conductivity fields, convec- 
tion fields, reaction fields, and so on. The dependence of the solution from these fields is nonlinear, while 
the solution linearly depends on external source terms / and boundary conditions. We assume there is a 
bilinear form associated with the operator L that allows us to write the variational form of ([1} in the form 

a(u,v;n) = (u),v), (3) 

for all test function v and where et(-,-;/i) is bilinear, coercive, and continuous for each /x, and (•,•) is an 
inner product. We assume that a(-, •; fi) is sufficiently smooth with respect to \x. 

The proposed algorithms have advantages when L^iu) has an affine representation defined as follows: 

Q 

L,(u):=J2^ML q (x). (4) 

1=1 

Here, L q {u) are heterogeneous operators with multiple scales with the coefficients that have high contrast, 
the parameter /i 6 A C R p is possibly a coarse-grid function and the functions <d q : A — > R. In terms of 
the corresponding bilinear form, we have 

Q 

a(u, v;n) = ^2 v). 

8=1 

This affine representation allows pre-computing coarse-scale projections of L q (x) in the offline stage and 
using them in the online stage. This reduces the computational cost considerably. 

Before discussing the GMsFEM, we introduce the notion of a coarse grid. Let T H be a usual conforming 
partition of D into finite elements (triangles, quadrilaterals, tetrahedrahals,...). We call this partition the 
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coarse grid and assume that this coarse grid is partitioned into fine-grid blocks. Each coarse-grid block is a 
connected union of fine-grid blocks. We denote by N v the number of coarse nodes, by the vertices 

of the coarse mesh T H and define the neighborhood of the node Xi by 

Wi = \J{Kj G T H ] XiGKj} (5) 

(see Figure 01 . 

The GMsFEM has a structure similar to that of MsFEM. The main difference between the two ap- 
proaches is that we systematically enrich coarse spaces in GMsFEM and generalize it by considering an 
input space consisting of parameters and source terms. In the first step of GMsFEM (offline stage), we 
construct the space of 'snapshots', hots , a large dimensional space of local solutions. In the next 

step of the offline computation, we reduce the space hots via some spectral procedure to V£g ■ In the 

second stage (online stage), for each input parameter, we construct a corresponding local space, V£j* that 
is used to solve the problem at the online stage for the given input parameter. Our systematic approach 
allows us to increase the dimension of the coarse space and achieve a convergence. 



2.1 An illustrative example 

Before presenting details of GMsFEM, we will present a simple example demonstrating the main concept 
of GMsFEM. We consider 

— (Hv(k(x; u)Vu) = f in D, 

u = on 3D and assume uq < u(x) < un, where uo and un are pre-defined constants. We assume the 
interval [uq,un] is divided into N equal regions uo < u\ < ... < ujv— l < un- 
As for the space of snapshots, Snapshots' we can consider 

— div(«(a;; Mj)V-0™ ap ) = in u>i, 

V , z"J- ap = ^i(x) on duji. Here, Si(x) are some set of functions defined on dui, e.g., unit source terms. We can 
also consider the space of fine-grid functions within tOj as the space of snapshots. 

As for offline space, V£g , we perform a spectral decomposition of the space of snapshots. We re- 
numerate the snapshot functions in uii by ■)/; i snap . We consider 

(S oS :=)s mn = [ ^ K (z; U ;)V« aP -V< nap , {A off :=)(W = f ^Mx- ^T" , 

where t\ are some weights. Here, S oil = (s mn ), and A oS := (a mn ). The choices for the matrix A oS will 
be discussed later. To generate the offline space, V£g, we choose the largest M g eigenvalues (see later 
discussions on M a s) of 

zioff.T.ofF \ off Qoff vryoff 

in. rn rn 

and find the corresponding eigenvectors in the space of Snapshots by multiplication, Y^j ^°f^j nav , where 
are coordinates of the vector ^° s . Note that the offline space is computed across all u ,..., u^. 
To compute the online space for a given u q (the value around which the global problem is linearized), 
we consider an eigenvalue problem in V£g ■ Lets denote the basis of V£g by ip^ . We consider a spectral 
decomposition of V£g via 

S on = (s mn ) - / K (x;u q )Vi;°* ■ VVf , A° n = (a mn ) = [ K(x;u q )i>°X- 
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To generate the online space, V£g , we choose the largest M on eigenvalues of 

4on^on _ xon^on 

m m m 

and find the corresponding eigenvectors in the space of by multiplication, J2j ^ij'tj 8 ^ where \E'°j 1 are 
coordinates of the vector ^>° n , denote these basis functions by ?/>™. 

At the final stage, these basis functions will be coupled via a global formulation, e.g., Galerkin for- 
mulation. In this case, the eigenfunctions are multiplied by the partition of unity functions to obtain a 
conforming basis. In this nonlinear example, one can use an iterative Picard iteration at the previous value 
of the solution u n (x), and in each iteration, a global problem is solved with V on for the value of u n (x) 
averaged over a coarse block. 

We will discuss the use of affine forms of k(x; u) and other details next. Also, these basis functions can 
be multiplied by a partition of unity or used in discontinuous Galerkin framework. Furthermore, various 
choices of spectral decomposition approaches can be employed. These details are also provided next. 



2.2 Step 1. Local multiscale basis functions (offline stage) 

1.1. Generating snapshots 

We call the local spatial fields (or local solutions ) as "local snapshots" . These local spatial fields are 
used to construct the offline space. This space consists of spatial fields defined on a fine grid, i.e., they are 
vectors of the dimension of fine-grid resolution of the coarse region. A common option for local snapshots 
is to use a fine-grid space (e.g., fine-grid linear functions) that resolves the coarse-grid block. However, in 
some cases, one can consider smaller and more appropriate spaces to construct the space of local snapshots. 
This will be discussed next. 

In the first step, a local reduced-order approximation is constructed based on the input space. Here, we 
will discuss two approaches and emphasize only the first approach where we will construct local approximate 
solutions by assuming that source term / is a smooth function. In this construction, the source term will 
not enter in the space of local solutions and its effect will be captured at the coarse-grid level via a global 
coupling. 

Remark 1. One can often ignore lower order terms in the operator L (see an example below) and use an 
operator different from the original one on a coarse grid to construct snapshots. To formalize this step, we 
assume that there exists such that if 

L^(u) — and u — u on dK 

then 

\\u-u\\ < 5 H , 

where Sn^OasH^O. In this case, we have a corresponding bilinear form a(u,v;/j,). For example, if 
L(u) = — div(n(x/ e, n)Vu + q(x, n) ■ Vu, one can use L(u) — —div(K(x/e,/j,)'Vu (see e.g., [Hi), provided, 
e.g., q is a bounded function. To avoid a cumbersome notation, we do not use L in the rest of the 
presentation. 

We construct local snapshots of the solutions that approximate the space of solutions generated by 

LM = (6) 
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in each subdomain Wi (see Figure [4]). We will consider two choices, though one can consider other options. 
First choice. We consider snapshots generated by 



= in w 4 (7) 

with boundary conditions 

$3 = b t in du)i, (8) 

with 6; being various source terms along the external boundary of u>i . One can also use Neumann boundary 
conditions or boundary conditions 6/ defined on larger domains for generating snapshots. 

Second choice. We can use local fine-scale spaces consisting of fine-grid basis functions within a 
coarse region. In this case, the offline spaces will be computed on a fine grid directly. The local fine-grid 
space has an advantage if the dimension of the local fine spaces is comparable to the dimension of K^ apshots 
computed by solving local problems as in the first choice. 

One can also construct local snapshots by solving 

L nMlj) = 9i in Ui 

with ipfi = on duji and where gi is supported on dK. The space spanned by ipf* is V^"* hots (cf. [T5]). 

Remark 2. On generating snapshots via the approximation of source term /. The snapshot 
spaces generated above can be larger than the space corresponding to local snapshots obtained from 

= / 

where f runs over the input space corresponding to the source term. Here, we take the restriction of Vf 
in uji to generate the space of local snapshots. This space can be taken as a span of (f)f, where (f)f are 
restrictions of the solutions of 

LM) = 4 

onto a coarse region uii . Here (f>{ are local basis functions that represent 

i 

These basis functions are global; however, they can be localized at a coarse-grid level (see e.g., Owhadi and 
Zhang, 2011, 14 ^j ) and we need a further dimension reduction of the space following Step 1.2. One can 
use the restriction of these global snapshots on the boundary of the coarse-grid block to compute the space 
of snapshots. This space of snapshots is an appropriate space of functions for which spectral decomposition 
needs to be performed. 

1.2. Generating the offline space 

Once the space of snapshots, hots = Span(^pf'-), is constructed for each u>i, spectral approaches 

are needed to orthogonalize and possibly reduce the dimension of this space. As a result, we will obtain 
the offline space, V^J that will be used to construct multiscale basis functions in the online stage. 

To perform a dimension reduction, we consider an auxiliary spectral decomposition of the space 
Snapshots = ^P an (' t Pi',j) ^ 0T eac ^ u i- Orrr objective is to construct a possibly small dimensional space 
V£g and use it for constructing multiscale basis functions for each /i in the online stage. In general, we will 
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seek the subspace V"^ such that for any fj, and ip G Snapshots (/-0 ( Snapshots 0-0 * s tne s P ace 01 snapshots 
which are computed for a given /i) , there exists ipo G SJg , such that 

a°^-Vo,^-V'o;M)dK!(^-V'o,V'-^o;M) ) (9) 

where a°*(</>, </>; //) and s°f (</>, (f>; /x) are auxiliary bilinear forms. In computations, this involves solving an 
eigenvalue problem with a mass matrix and the basis functions are selected based on dominant eigenvalues. 
Note that this eigenvalue problem is formed in the space of snapshots. We propose two procedures for 
constructing V^g . We will discuss the requirements for the space Ss later on. 

Remark 3. In general, affl and contain partition of unity functions, penalty terms, and other dis- 
cretization factors that appear in coarse-grid finite element formulations. The norm corresponding to 
needs to be stronger, in general, to allow the decay of eigenvalues (see Remark^. However, one can also 
take s°$ to be weaker. 

We consider two other options. 

Option 1. We consider a°*(0, <fi; /z) and s^.{4>, 4>) f) to be independent of /x, i.e., 
afftf, cj>; /i) = <(</>, 0) and S ^(0, 0; = s$(<f>, 4>) . 

In this case, finding reduces to performing spectral decomposition of Snapshots with corresponding- 
inner products. A subspace V"^ of Snapshots such that for each ip 6 Snapshots' there exists ip Q G Ss such 
that 

a° ff (V -M- *M r< Ss$(il> -M- *M (10) 

for some prescribed error tolerance 5. We give some examples for the elliptic equation @. 
Example 1. In this example, we average the parameter /i to obtain a oS and s off . 



3 3 Jui 



(11) 



where t» are non-negative weights (a case with a fixed value of fj, is a special case) and Xi is a partition of 
unity function supported in w^. 

Example 2. For any ip G Snapshots' ^ = Y,i, 3 a i,j^i'j ( see ©)' we define ^ = Y,i a i,j^i'j- Not e that 
V-> = X) j • We use 

, (12) 
<W)=E^ / «(^Mj)|V(Xi^-| 2 , 

where tj are non-negative weights. 

To formulate the eigenvalue problem corresponding to (|TU]) . we will re- numerate the basis, Snapshots = 
Span(ijjj'). Then (|10[) will yield an algebraic eigenvalue problem 
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Here, A oS = (A° S j) and S oS = (S°j) are corresponding local matrices that are given by 

Afj = afS^T^?), Sfj = sZ^T^T), 

ip'j'itpj* G Snapshots- Here, we assume that is a non-degenerate positive definite matrix. Note that the 
matrices A oS and S° s are computed in the space of snapshots. The space V£g is constructed by selecting 
Li eigenvectors corresponding to largest eigenvalues. If we assume that 

Af > ... > Xf, 

then we choose the first Li eigenvectors in each uii corresponding to the largest Li eigenvalues to form the 
offline space. In particular, the corresponding eigenvectors are computed by multiplication, J^j ^if^ 1 m 
each u>i, where are coordinates of the vector \&° s . If we use the second example, then the matrices 
A oS and S oS are block diagonal with one block for each fij . 

Option 2. The second option for a numerical setup of ((5J) is the following (see [29]). First, we will 
compute the local spaces for each /i, V^h (jjl) such that for any /i and if) G Snapshots 0") (^ e -' snapshots 
are computed for a given /it), there exists ipo G (//), such that 

afiyi - Vo, V> - V'o; m) d -M-fhw)- (13) 

Secondly, we seek a small dimensional space V^g (across all /i's) such that any t/^o G V"}* „(/i) can be 
approximated by an element of V^j. One can do it by defining a distance function between two spaces 
Vgff and V"^ ) „(a*2) and identifying a small number of /i's such that the spaces V^hjji) corresponding 
to these /i's approximate the space spanned by (/i) for all /i's. Then, the elements of these V ^ M (/i)'s 

will form the space . 

Note that it is important to have a small dimensional space of snapshots to reduce the computational 
cost that arises in the calculations of multiscale basis functions. The result of Step 2 is the local space V£g. 

Remark 4 (On S norm). In this remark, we show eigenvalues decay for several choices of and . 
We consider a target domain to t = [0.4,0.6] x [0.4,0.6] within D — [0,1] x [0,1]. Our space V snaps hots is 
generated by solving problems with force terms located outside w t - We consider several choices for and 
s°ff . First, we choose 



. 1 Jolt Jul* 



where Xi o,re bilinear basis functions. We also consider 

■,2,off_ 



a ui t 



k(x)\u\ 2 , slf= \ k(x)\Vu 



We also consider 



a 3 jf= / k(x)\Vu\ 2 , 4;f = / K (x)\Vu\ 2 , 



where uj ext = [0.3,0.7] x [0.3,0.7]. The last choice is motivated by Babuska and Lipton, 2011, JlOf). where 
a larger domain uj^ , oj{ C oj^~ (see Figure^, is used for eigenvalue computations. 
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Figure 3: Top left: Force locations. Top right: Eigenvalues when using a ' and s > . Bottom left 
Eigenvalues when using a 2 ' off and s 2 ' off . Bottom right: Eigenvalues when using a 3,oS and s 3 ' off . 

We take k(x) to be 10 2 in [0.45,0.55] x [0.45,0.55] and 1 elsewhere. In Figure^ we plot the eigenvalues 
corresponding to the above choices of and s°^ . The space of snapshots are generated using unit force 
terms in the locations shown in Figure^ top left. As we see in all the cases the eigenvalues decay and one 
can choose a threshold for the cut-off. The decay of eigenvalues depends on the choices of and , 
and also on the choice of \i ■ 

We note that one can combine the snapshot calculations with a dimension reduction by using greedy 
algorithms instead of solving an eigenvalue problem as discussed above. This will provide an estimate in 
sup-norm with respect to \x. To illustrate this we consider the offline quotient of Example 2 presented 
before. In this example we need to 

• select adequate values among all possible values of the parameter \i 

• compute the important modes of this eigenvalue problem in Example 2. 

We can use the Reduced Basis (RB) methodology to achieve these objectives iteratively. This is due to 
the fact that the eigenvalue computation for each \Xj in Example 2 is decoupled from each other. The 
whole set of eigenvalue problem can be seen as a parameter dependent problem and a RB methodology 
can be used. More precisely, we would like to choose a set of parameters {fi m } such that it minimizes the 
error between the true important modes and their approximations in some metric |j • \\ as associated with 
s°^(-0, -0)- We seek a set of {yu m } such that 




inf sup 1 1 ipf -<pl 

{Mm} \fieA 



II 



as 
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where ip^ % ' app is an approximation to the true important mode tp^*. In practice, we do not compute 
i^ i {ix),Li G A for two reasons: first A can be a continuous set, second computing ^'(/x) for a given fx is 
computationally expensive, and so we just want to perform this computation for the selected {/i m }. The 
RB approach suggests a feasible procedure, which consists of selecting the minimizers of 

{fi m } = inf sup Ajv(aO ) j 

(C'"} W Atrial / 

where Atrial is a discrete set and Ajy(/x) is an estimator of the quantity e(/i). We refer to |29j for details 
regarding Example 2. 

In conclusion of Step 1, we obtain the space of offline multiscale basis functions, V"g. The requirement 
is that the dimension of V^g is less (or equal to) than the number of fine-grid blocks within Wj. 

2.3 Step 2. Computing online multiscale basis functions and their coupling 

At the online stage, for each parameter value, multiscale basis functions are computed based on each local 
coarse region. In particular, for each and for each input parameter, we will formulate a quotient for 
finding a subspace of (fx) where the space will be constructed for each ix (independent of source terms). 
We seek a subspace V£* (li) of such that for each tp € V£g , there exists ipo € V££ (fx) such that 

a ^-4^-*;rf^CW-*i-*;/') (14) 

for some prescribed error tolerance 6 (different from the one in the offline stage), and the choices of a° n 
and s° n . The corresponding eigenvalue problem is formed in the space of offline basis functions. We note 
that an assumption as in Remark [3] (see also Remark [J) is needed for obtaining a convergence result and, 
in general, a° n and s° n contain partition of unity functions, penalty terms, and other discretization factors 
that appear in finite element formulations. In particular, we assume that ok {Xi^i Xi^', A*) ^ a ™("07 "0; m) 
(where <xk corresponds to ^ in K). 

Remark 5. The bilinear form a° n can be chosen based on finite element variational formulation in u>i (see 
Remark^. In particular, it can be chosen to be the same as the coarse-grid finite element formulation 
(see U8\) . U9\) , and (B(J\) ). As before, one can also take s°-^ to be weaker, provided the error in s-norm 
between the solution and its approximation can be estimated. 

We note that a choice of partition of unity is important to achieve smaller dimensional coarse spaces 
and one can look for a choice of optimal partition of unity function. Because is a finite dimensional 
space, the following local eigenvalue problem 

A° n vr = \f a s on i/)? n (is) 

can be used to find the basis functions, where A on and 5 on are local matrices corresponding to a(-, •; fx) 
and s(-, •; fx), A on = (A™) and S on = (S°j) are corresponding local matrices that are given by 

Here, we assume that S 1 ™ is a positive definite matrix. Note that the matrices A on and S on are computed 
in the space of offline basis functions. (fx) is constructed by selecting Li eigenvectors corresponding 
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to the largest eigenvalues. In particular, the corresponding eigenvectors are computed by multiplication, 
J2j ^i^j 1 m eacn w i' where x P° n are coordinates of the vector \P° n . 

We need a certain requirement for the space to guarantee that (fTTJ|) holds for all tp £ Snapshots ( anc ^ 
not only for all ip £ which may be a smaller subspace). Indeed, we need to guarantee that the space 
Km wm have approximation properties similar to those of the space V££ computed using all snapshots 
Snapshots (instead of V£g that may be smaller subspace). In particular, we want to guarantee that (jT5J) 
holds for any if) € Snapshots w ^ n ^0 G Kn* • One can show that under the conditions 

s$W, i/O i CWs 1>\ m). W> G CWs (16) 

the inequality (flU)) holds for all ip G V^"* hots . This can be shown using the triangle inequalities and (ITC1) . 

Once multiscale basis functions are constructed, we project the global solution onto the space of basis 
functions. One can choose different global coupling methods and we present some of them. 

Galerkin coupling. For a Galerkin formulation, we need conforming basis functions. We modify V££ 
by multiplying the functions from this space with partition of unity functions. The modified space has the 
same dimension and is given by Spanj{xii)^ l '° n ), where ip^ z ' on e V^CaO an d Xi is supported in Wj. Then, 
the Galerkin approximation can be written as 



u G 

If we introduce 



V G = Spamjixi^n, (17) 



then Galerkin formulation is given by 

a(u« 1 ,»; / i) = (/,») ) V«et (18) 

Petrov-Galerkin coupling. We denote V PG = Sparii^^^} and write the PG approximation of the 
solution as 

i,3 



a ms ' 



Then the Petrov-Galerkin formulation is given by 

a(u PG ,v;ti) = (/,«), V^t (19) 

where Vg is defined with ([T7|) . 

Discontinuous Galerkin coupling 

One can also use the discontinuous Galerkin (DG) approach (see also [HI [331 [5U] ) to couple multiscale 
basis functions. This may avoid the use of the partition of unity functions; however, a global formulation 
needs to be chosen carefully. We have been investigating the use of DG coupling and the detailed results 
will be presented elsewhere. Here, we would like to briefly mention a general global coupling that can be 
used. The global formulation is given by 

a DG (u,v) = f(v) for all v = {v K eV£], (20) 
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Figure 4: Schematic description of oversampled regions. 



where 



a DG {u,v)=Y, a K G M and /W = E/ f v « dx ( 21 ) 

K K ^ K 



for all u = {uk}, v = {vk}- Each local bilinear form a^ G is given as a sum of three bilinear forms: 

a K G (u, v) := a K (u,v) + r K (u,v) + p K (u,v), (22) 

where ok is the bilinear form, 

clk(u, v) := / kkVuk-VvkcLx, (23) 
Jk 

where kk is the restriction of k(x) in K; the tk is the symmetric bilinear form, 



r K (u,v) := j- [ Ke ( ^^-{vk -v K ') + ^^(uk' -Uk)) ds, 

e ^ k 1 eJe \dn K dn K ) 



where is a weighted average of k(x) near the edge E, Ie is the length of the edge E, and K' and K are 
two coarse-grid elements sharing the common edge E; and px is the penalty bilinear form, 

p K (u,v):= V] t^e I k e (u k > - u K )(v K / - v K )ds. (24) 

Here 5b is a positive penalty parameter that needs to be selected and its choice affects the performance of 
GMsFEM. One can choose eigenvalue problems based on DG bilinear forms. 

Remark 6. If the norm of interest is other than a(-, ■; p), we can use the Petrov-Galerkin approaches with 
a particular choice of test functions (e.g., \22j). 

Remark 7 (On the convergence of GMsFEM). Here, we briefly discuss how spectral convergence can enter 
in the error analysis for GMsFEM based on the Galerkin coupling. We denote the bilinear form f3J) on the 
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whole domain D by an(u,v; (j,), and u ms = J2i j c i.jXi4>) , an d u o % = 2j c t,j'V'}- Then, the error analysis 
is the following. 

a D (u - u mS) u- u ms ;n) = a D (^2xi(u - Up*) Xi(u ~ u o');^) r< 

i i 
i i i 

i i i 

ffere, we used f/ie /act i/iai axixi^^ XiV'i A 4 ) ^ ^"(XiVs XiV'j A*); which is an assumption on the choice ofa° n , 
the inequality H0\) , and replaced the sum over K by the sum over larger regions Wj. Moreover, we have used 
the fact that the solution restricted to Ui is in the space of offline basis functions. As we mentioned, one can 
choose a° n (xi'0, Xi"0; A 4 ) based on finite element formulation such that a°"(xiV': Xi^l A 4 ) = a w, {Xi^i Xi^'i A*)- 
T7ie hidden constants depend on the number of neighboring element of each coarse block and the constants 
in A10\) . |A^| ) and M6\) . We note that one can obtain sharper estimates using a bootstrap argument in some 
cases (see J30\j) and show that the error decreases as the coarse-mesh size decreases. 

2.4 Iterative solvers - online correction of fine-grid solution 

In the previous approach, the coarse spaces are designed to achieve a desired accuracy. One can also 
iterate (on residual and/or the basis functions) for a given source term and converge to the true solution 
without increasing the dimension of the coarse space. Both approaches have their application fields and 
allows computing the fine-grid solution with an increasing accuracy. Next, we briefly describe the solution 
procedure based on the concept of two-level iterative methods. 

We denote by {D' i }f^ 1 the overlapping decomposition obtained from the original non-overlapping de- 
composition {Di}^^ by enlarging each subdomain Di to 

D'i = Di U {x e £>,dist(x, A) < Si}, i=l,...,M, (26) 

where dist is some distance function and let V^ l (D' i ) be the set of finite element functions with support 
in D\ and zero trace on the boundary dD[. We also denote by Rf : V" ' l (Z?^) — > V h the extension by zero 
operator. 

With the help of (u — uq), we can correct the coarse-grid solution. A number of approaches can be 
used. For instance, we write the solution in the form 

u = u Q + ^2xiVi, 

i 

where v\ are defined in u>i (though it can be taken to be supported in a different domain) and Xi i s a 
partition of unity. Suppose that Uj has zero trace on du>i. We can solve the local problems 

L ( v i) = f - L fi(uo), 

with zero boundary condition. Other correction schemes can be implemented to correct the coarse-grid 
solution. For example, we can use the traces of uq to correct the solution in each subdomain in a consecutive 
fashion. 
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When the bilinear form is symmetric and positive definite, we can consider two-level domain decompo- 
sition additive methods to find the solution of the fine-grid finite element problem 

a(u, v; fj) = f(v), for all v £ V h , (27) 

where V is the fine-grid finite element space of piecewise linear polynomials. The matrix of this linear 
systems is written as 

A(//>(» = b. 

Here, A is the stiffness matrix associated to the bilinear form a and b such that v T b = f(v) for all v £ V ■ 
We can solve the fine-scale linear system iteratively with the preconditioned conjugate gradient (PCG) 
method (or any other Krylov type method for a non-positive problem). Any other suitable iterative 
scheme can be used as well. We introduce the two level additive preconditioner of the form 

M 

B-\») = Rl on {^)A^(p)R a .on(y) + J2 R I A i 1 (») R i> (28) 

i=i 

where the local matrices are defined by 

v T Ai(n)w = a{v,w,(j,) for all v,w £ V£(D'A. (29) 

The coarse projection matrix i?^ on is defined by R^ on — Ro on (^) — [$i, ■ ■ • , &N V ] an d the online coarse 
matrix A (fi) — i?o i0n (/x).A(/i)i?o on (/x). The columns <IVs are fine-grid coordinate vectors correspond- 
ing to the basis functions, e.g., in the Galerkin formulations they correspond to the basis functions 
{Xi{ x )i )l j % ' 0n { x 'i A 1 )}- See [541 146] and references therein for more details on various domain decomposi- 
tion methods. 

The application of the preconditioner involves solving local problems in each iteration. In domain 
decomposition methods, our main goal is to reduce the number of iterations in the iterative procedure. It 
is well known that a coarse solve needs to be added to the one level preconditioner in order to construct 
robust methods. The appropriate construction of the coarse space Vb plays a key role in obtaining robust 
iterative domain decomposition methods. Our methods provide an inexpensive coarse solves and efficient 
iterative solvers for general parameter-dependent problems. This will be discussed in the next sections. 

3 Case studies and relation to existing methods. Discussions and 
applications 

In this section, we illustrate basic concepts via some specific examples. We use existing methods in the 
literature for multiscalc problems and show how these methods can be put under the general framework 
of the GMsFEM and show some numerical results. 

3.1 Case with no parameter 

In this section, we write a method proposed in [25] as a special case of GMsFEM. First, we consider a case 
with no parameter, i.e., 

L(u) = /, or corresponding a(u, v) — f(v). 
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In this case, offline multiscale basis functions are used for online simulations due to the absence of the 
parameter. Next, we discuss the construction of multiscale basis functions (see [251127] ). 

The construction of the offline space starts with the space of snapshots. We choose the space of 
snapshots as all fine-grid functions within a coarse region. 

For the construction of the offline space, we can choose 

<W)= E a "dXk*P,Xki>), (30) 

k,u> k (~)u>iy£0 

where a Wi (-, ■) is the restriction of a(-, ■) in Uj and Xk is the partition of unity corresponding to lo^. For 
example, for the elliptic equation (see (J2j), 

As for s°*, we select 

The corresponding eigenvalue problem can be explicitly written as before. In our numerical implementation, 
we choose a°*(^, ip) — J2k u k n^#o JL K I^Xfe| 2 l' ( /'| 2 ( see [25] that shows that this is not smaller than ([50)1 
in the space of local k- harmonic functions) . 

In this example, we do not construct an online space and the offline space is used in GMsFEM, 
Km = VoS ■ The basis functions are constructed by multiplying the eigenvectors corresponding to the 
dominant eigenvalues by partition of unity functions, see ([P7|) . The Galerkin coupling of these basis 
functions are performed based on (|18[) . Note that the stiffness is pre-computed in the offline stage and 
there is no need for any stiffness matrix computation in the online stage. 

We point out that the choice of initial partition of unity basis functions Xi > are important in reducing 
the number of very large eigenvalues. We note that the dimension of the coarse space depends on the 
choice of Xi an d, thus, it is important to have a good choice of Xi- The essential ingredient in designing 
them is to guarantee that there are fewer large eigenvalues, and thus the coarse space dimension is small. 
With an initial choice of multiscale basis functions Xi that contain many localizable small-scale features of 
the solution, one can reduce the dimension of the resulting coarse space. 

Next, we briefly discuss a few numerical examples (see [25] for more discussions). We present a numerical 
result for the coarse-scale approximation and for the two-level additive preconditioner ([28]) with the local 
spectral multiscale coarse spaces as discussed above. The equation — div(KVw) = 1 is solved with boundary 
conditions u = x + y on dD. For the coarse-scale approximation, we vary the dimension of the coarse 
spaces by adding additional basis functions corresponding to the largest eigenvalues. We investigate the 
convergence rate, while for preconditioning results, we will investigate the behavior of the condition number 
as we increase the contrast for various choices of coarse spaces. The domain D = [0, 1] x [0, 1] is divided 
into 10 x 10 equal square subdomains. Inside each subdomain we use a fine-scale triangulation, where 
triangular elements constructed from 10 x 10 squares are used. We consider the scalar coefficient k(x) 
depicted in Figure [5J that corresponds to a background one and high conductivity channels and inclusions. 

We test the accuracy of GMsFEMs when coarse spaces include eigenvectors corresponding to the large 
eigenvalues. We implement GMsFEM by choosing the initial partition of unity functions to consist of 
multiscale functions with linear boundary conditions (MS) (see ([3Tjl . We use the following notation. 
GMsFEM+0 refers to the GMsFEM where the coarse space includes all eigenvectors that correspond to 
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H = 1/10 




I- Ik 


\-\h 


A* 


GMsFEM+0 


14.3 


14.3 


6.17e-02 


0.0704 


GMsFEM+1 


5.82 


5.82 


1.02e-02 


0.0117 


GMsFEM+2 


5.33 


5.33 


8.60e-03 


0.0071 


GMsFEM+3 


4.64 


4.64 


6.53e-03 


0.0043 


GMsFEM+4 


4.01 


4.01 


4.80e-03 


0.0032 



Table 1: Convergence results (in %) for GMsFEM with MS with increasing dimension of the coarse space. 
Here, r\ = 10 6 . The initial coarse space is spanned by multiscale basis functions with piecewise linear boundary 
conditions (x ms )- The coefficient is depicted in Figure [5j 



eigenvalues which are asymptotically unbounded as the contrast increases, i.e., these eigenvalues increase as 
we increase the contrast. One of these eigenvectors corresponds to a constant function in the coarse block. 
GMsFEM+n refers to the GMsFEM where in addition to eigenvectors that correspond to asymptotically 
unbounded eigenvalues, we also add n eigenvectors corresponding to the next n eigenvalues. 

In previous studies [34] [35j [30], we discussed how the number of these asymptotically unbounded 
eigenvalues depends on the number of inclusions and channels. In particular, we showed that if there are 
n inclusions (isolated regions with high conductivity) and m channels (isolated high-conductivity regions 
connecting boundaries of a coarse grid), then the number of asymptotically unbounded eigenvalues is 7i + m 
when standard bilinear partition of unity function, %?, used. However, if the partition of unity Xk is chosen 
as multiscale finite element basis functions ([38]) defined by 

div( K VxD = in K e Ui, xT = Xi in dK, V K e ua, (31) 

then the number of asymptotically unbounded eigenvalues is m. We can also use energy minimizing basis 
functions that are defined (see [57]) as 

min£ f n\VxT S \ 2 (32) 

subject to Yli Xi™^ = 1 with Supp(xj) C u>i, i = 1,...,N V , to achieve even smaller dimensional coarse 
spaces. 

In all numerical results, the errors are measured in the energy norm (| • |^), H 1 norm ( | • |^i), and 
L 2 -weighted norm (| • \ 2 L 2 ) respectively. We present the convergence as we increase the number of additional 
eigenvectors. In Table [T] we present the numerical results when the initial partition of unity consists of 
multiscale basis functions with linear boundary conditions for the contrast r\ = 10 6 . We note that the 
convergence is robust with respect to the contrast and the error reduces. The error is proportional to the 
largest eigenvalue (A*) whose eigenvector is not included in the coarse space as one can observe from the 
table (correlation coefficient between A* and the energy error is 0.99). We observe that the errors are 
smaller compared to those obtained using MsFEM with piecewise linear initial conditions. 

Next, we present the results for a two-level preconditioner. We implement a two-level additive pre- 
conditioner with the following coarse spaces: multiscale functions with linear boundary conditions (MS); 
energy minimizing functions (EMF); spectral coarse spaces using piecewise linear partition of unity func- 
tions as an initial space (GMsFEM with Lin); spectral coarse spaces where multiscale finite element basis 
functions with linear boundary conditions (x ms ) ar e used as an initial partition of unity (GMsFEM with 
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MS); spectral coarse spaces with k where energy minimizing basis functions (x e ) are used as an initial 
partition of unity (GMsFEM with EMF). In Tabled we show the number of PCG iterations and estimated 
condition numbers. We also show the dimensions of the coarse spaces. Note that the standard coarse space 
with one basis per coarse node has the dimension 81 x 81. The smallest dimension can be achieved by 
using energy minimizing basis functions as an initial partition of unity. We observe that the number of 
iterations does not change as the contrast increases when spectral coarse spaces are used. This indicates 
that the preconditioncr is optimal. On the other hand, when using multiscale basis functions (one basis 
per coarse node), the condition number of the preconditioned matrix increases as the contrast increases. 




Figure 5: Coefficient k(x). The dark region has high conductivity r\ and the white background has conduc- 
tivity 1. Green lines show the coarse grid. 



V 


MS 


EMF 


GMsFEM with Lin 


GMsFEM with MS 


GMsFEM with EMF 


10 3 


83(2.71e+002) 


69(1.43e+002) 


31(8.60e+000) 


31(9.34e+000) 


32(9.78e+000) 


10 5 


130(2.65e+004) 


74(1.29e+004) 


33(8.85e+000) 


33(9.72e+000) 


34(1.02e+001) 


10 7 


189(2.65e+006) 


109(1.29e+006) 


34(8.85e+000) 


35(9.60e+000) 


37(1.02e+001) 


Dim 


81 


81 


165 


113 


113 



Table 2: Number of iterations until convergence and estimated condition number for the PCG and different 
values of the contrast r\ with the coefficient depicted in Figure We set the tolerance to le-10. Here 
H = 1/10 with h = 1/100. 

3.2 Elliptic equation with input parameter 

Here, we present a method proposed in 29 as a special case of GMsFEM. We consider a parameter- 
dependent elliptic equation (see ©) 

L^(u) = /, or corresponding a{u,v] fi) = f(v). (33) 

As before, we start the computation with the space of snapshots consisting of local fine-grid functions 
and compute offline basis functions. In this example, we will construct offline multiscale spaces for some 
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selected values of /i, /ij (i = 1, -/V r {,), where N r b is the number of selected values of \x used in constructing 
multiscale basis functions. Note that these values are chosen automatically using RB technique. These 
values of \i are selected via an inexpensive RB procedure [25]. We briefly describe the offline space 
construction. For each selected fj,j, we choose the offline space as in the previous example, 

where a UJi {ip, ip; n) is the restriction of a(ip, ip; /i) onto and \k is the partition of unity corresponding to 
u>k- For example, for elliptic equation, 



J Di 

As for s°*, we select 

Then, the selected eigenvectors are orthogonalized with respect to H 1 inner product. In our numerical 
implementation, we choose a^j(tp, 0) = 2~Zfc,u fc n^#o Xj< k ( x > 

The resulting a^(ip, ip) and s^(ip, ip) can also be defined in the following way. For any ip G V^"* hots' 
V-> = J2i j a i,j^?p we define ipj = ctijipfj. Note that -0 = ^ Vj- Then, s oS and a off are defined in the 
following way. 

where are non-negative weights. In this case, the weights do not play any role as the resulting algebraic 
eigenvalue problem is block diagonal. 

At the online stage, for each parameter value, multiscale basis functions are computed based on the 
solution of a local problem. In particular, for each cjj and for each input parameter, we formulate a quotient 
to find a subspace of V^(fi), where the space will be constructed for each For the construction of the 
online space, we choose 

k,io k (~) tdi/O 

Here \k is the partition of unity corresponding to uik and a UJi (ip, ip; /i) is the restriction of a(tp, ip] fi) onto 
u>i. For the bilinear form for s, we choose 

In this case, the online space is a subspace of the offline space and computed by solving an eigenvalue 
problem for a given value of the parameter /i. The online space is computed by solving an eigenvalue 
problem in w.; using V^f, see (fT^j) . Using dominant eigenvectors, we form a coarse space as in (fT?)) and 
solve the global coupled system following (fT5)l . For the numerical example, we will consider the coefficient 
that has an affine representation (see The online computational cost of assembling the stiffness matrix 
involves summing Q pre-computed matrices corresponding to coarse-grid systems. 
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We point out that the choice of initial partition of unity functions, Xii are important in reducing the 
number of very large eigenvalues (we refer to [29. for further discussions). 

We present a numerical example for — 6iv{k{x] /i) Vu) = 1 which is solved with boundary conditions 
u = x + y on 3D. We take D = [0, 1] x [0, 1] that is divided into 10 x 10 equal square subdomains. As 
in Section 13.11 in each subdomain we use a fine-scale triangulation, where triangular elements constructed 
from 10 x 10 squares are used. 

3.2.1 One dimensional input parameter 

We consider a permeability field which is the sum of two permeability fields with each containing inclusions 
such that their sum gives a channelized permeability field. The permeability field is described by 



k(x; fi) := (1 - ^i)ko(x) + ^lk 1 (x). 



(34) 



We can represent 3 distinct different features in n(x; fi): inclusions (left), channels (middle), and shifted 
inclusions (right), see Figure |51 There exists no single value of /x that has all the features. Furthermore, 
we will use a trial set for the reduced basis algorithm that does not include fi = 0.5. 



Figure 6: From left to right: \i = 0, fj, = 1 and /i = 1/2. 

As there are three distinct spatial fields in the space of conductivities, we will choose several functions 
in our reduced basis. In the case of insufficient number of samples in the offline space, we will observe 
that when the online permeability field does not contain appropriate features, then we can not obtain the 
convergence. We observe in Table |3] that we indeed need N r b > 3 to capture all the details of the solution. 
In these tables, we compare the errors obtained with GMsFEM when the online problem is solved with a 
corresponding number of basis functions. We observe a convergence with respect to the number of local 
eigenvectors when N r b is chosen such that it contains spatial features included in the conductivity space. 
The error is proportional to the largest eigenvalue whose eigenvector is not included in the coarse space as 
one can observe from the table (correlation coefficient between A* and the energy error is 0.99). We have 
also computed weighted Li error which shows a similar trend and these errors are much lower. 

We use this conductivity field example to test the two-level additive Schwarz preconditioners. In this 
example, we only use coarse spaces based on reduced models. The numerical results are presented in 
Table [4] In these numerical results, we observe that when the dimension of the reduced space is 3 (and 
more), the condition number of the preconditioned system is independent of the contrast. In this example, 
we only choose the eigenvectors that correspond to asymptotically unbounded eigenvalues. Note that in 
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H = 1/10 


N rb = 2 


A* 


N rb =3 


A* 


GMsFEM+0 
GMsFEM+1 
GMsFEM+2 


65.20(121) 
18.87(462) 
15.37(566) 


0.4469 
0.0874 
0.0035 


25.92(235) 
7.22(356) 
6.07(477) 


0.0311 
0.0040 
0.0028 



Table 3: Convergence results (energy norm in % and space dimension) for GMsFEM with the increasing 
dimension of the coarse space. Here, h = 0.01, 77 = 10 6 , and y, = 1/2 (error with MsFEM 39.29%). 

this example, we only choose basis functions corresponding to the interior nodes, while in the coarse- 
grid approximation, we choose basis functions that also represent boundary nodes. We observe that the 
number of iterations does not change as the contrast increases when spectral coarse spaces are used. On 
the contrary, when using multiscale basis functions (one basis per coarse node), the condition number of 
the preconditioned matrix increases as the contrast increases. The latter is due to the fact that the coarse 
space does not contain enough degrees of freedom. 



V 


MS 


N rb = 2 


Nrb = 3 


10 4 


80(7.26e + 002) 


38(1. 76e + 001) 


33(9.64e + 000) 


10 6 


120(7.17e + 004) 


44(2. 64e + 001) 


35(9.85e + 000) 


Dim 


81 


310 


338 



Table 4: Number of iterations until convergence and estimated condition number for the PCG and different 
values of the contrast 77 and y = 1/2. We set the tolerance to le-10. Here H = 1/10 with h = 1/100. 

3.2.2 Four dimensional input parameter 

We consider a permeability field which is the sum of four permeability fields each of which contains 
inclusions such that their sum gives several channelized permeability field scenarios. The permeability 
field is described by 

k(x; h) := + ^iki(x) + ^3^3(2;) + haka,{x). (35) 

There are several distinct features in this family of conductivity fields which include inclusions and the 
channels that are obtained by choosing fj,i = fi2 = 1/2 or 7x3 — 714 — 1/2. There exists no single value of 
fi that has all the features. Furthermore, we will use a trial set for the reduced basis algorithm that does 
not include jiij = 1/2 (i = 1,2,3,4). 

As there are several distinct spatial fields in the space of conductivities, we will choose multiple functions 
in our reduced basis. In the case of insufficient number of samples in the offline space, we observe that the 
online permeability field does not contain appropriate features and this affects the convergence rate. We 
observe in Table [5] that we indeed need N r b > 4 to capture all the details of the solution. In this table, we 
compare the errors obtained by GMsFEM with a different number of online basis functions. We observe 
convergence with respect to the number of local eigenvectors when N r b increases. We note that in these 
computations, we also have an error associated with the fact that the offline space is not sufficiently large 
and thus the error decay is slow as we increase the number of basis functions. We have also computed 
weighted Li error which shows a similar trend and the Li errors are generally much smaller. 
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Figure 7: From left to right and top to bottom: Ki, K2, «3, K4. 



H = 1/10 


N rb = 2 


N rb = 3 


N rb =4 


GMsFEM+0 
GMsFEM+1 
GMsFEM+2 
GMsFEM+3 


45.5(478) 
39.4(599) 
38.4(720) 
36.2(841) 


40.3(565) 
27.5(686) 
26.8(807) 
26.2(928) 


9.2(588) 
5.3(709) 
5.1(830) 
4.9(951) 



Table 5: Convergence results (energy norm in % and space dimension) for GMsFEM with the increasing 
dimension of the coarse space. Here, h = 0.01, 77 = 10 6 , and A*i = A*2 = A*3 = A*4 = 1/2 (error with MsFEM 
96.77%). 



3.3 Anisotropic flows in parameter-dependent media 

In this example, we apply GMsFEM to anisotropic flows by considering the elliptic problem with tensor 
coefficients 

where the «i 1 coefficient is described by 

ki A (x; fx) := (1 - fx)K (x) + /i«i(ar). (36) 

We consider an example from Section [3 . 21 where the first component of the permeability has 3 distinct 
different features in k(x; /i): inclusions (left), channels (middle), and shifted inclusions (right); see Figure[6] 
In this example, there exists no single value of \i that has all the features. We will use a trial set for the 
reduced basis algorithm that does not include fi — 0.5. The trial set is chosen as in the case of isotropic case 
using a greedy algorithm. The important features in these permeability fields characterize the preferred 
directions of conductivity. 
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We note that the coarse space has a structure that differs from the case of isotropic coefficients. In 
particular, the coarse space has a larger dimension and contains all functions that are constant along x\ 
direction in the examples under consideration |28) . 

First, we present numerical results for the two-level domain decomposition solvers. In Table |H1 we 
present the results for a two-level preconditioner. We implement a two-level additive preconditioner with 
the coarse spaces introduced earlier: multiscale functions with linear boundary conditions (MS); spectral 
coarse spaces where piecewise linear partition of unity functions are used as an initial partition of unity 
(GMsFEM with Lin); and spectral coarse spaces where multiscale finite element basis functions with linear 
boundary conditions (x ms ) are used as an initial partition of unity (GMsFEM with MS). 

From this table, we see that the number of PCG iterations and estimated condition numbers do not 
depend on the contrast when spectral basis functions are used. We use N r b — 3 because we need at least 
three features to represent the permeability field as in the earlier example. We observed from our previous 
study that if N r b < 3, one could not get a contrast- independent condition number for the preconditioned 
system. On the other hand, when multiscale basis functions (one basis function per node) is used, the 
number of iterations and the condition number of the preconditioned system increase as the contrast 
increases. These results indicate that the preconditioner is optimal when the spectral basis functions are 
used and the coarse spaces include eigenvectors corresponding to important eigenvalues. Moreover, we 
observe that the dimension of the coarse space is smaller when multiscale finite element basis functions are 
used an initial partition of unity. 



V 


MS 


GMsFEM with Lin N rb = 3 


GMsFEM with MS N rb = 3 


10 4 


125(6.52e +2) 


45(22.79) 


42(15.16) 


10 6 


260(6. 14e +4) 


37(8.94) 


44(13.15) 


Dim 


81 (0.8% of fine DOF) 


862 (8.4% of fine DOF) 


744 (7.4% of fine DOF) 



Table 6: Number of iterations until convergence and estimated condition number for the PCG and different 
values of the contrast rj with fi = 0.5. We set the tolerance to le-10. Here H = 0.1 with h = 0.01. 

In Table [3 we present numerical results to study the errors of GMsFEM when the contrast is r] = 10 6 . 
As there are three distinct spatial fields in the space of conductivities, we choose at least three realizations. 
In Table El we compare the errors obtained by GMsFEM when the online problem is solved with a 
corresponding number of basis functions. We observe convergence with respect to the number of local 
eigenvectors. The convergence with respect to N r b can also be observed. 



H = 0.1 


N r b = 2 




N rb = 4 


GMsFEM+0 
GMsFEM+1 
GMsFEM+2 


3.8(1204) 
2.6(1325) 
2.1(1446) 


3.4(1360) 
2.4(1481) 
1.8(1602) 


3.2(1517) 
2.2(1638) ) 
1.7(1759) 



Table 7: Convergence results (energy norm in % and space dimension) with the increasing dimension of the 
coarse space. Linear basis functions are used to generate the mass matrix. Here, h = 0.01, rj = 10 8 , \i = 1/2 
(error with MsFEM 88%). 
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3.4 Parabolic problems 

Next, we discuss the extension of our procedures to parabolic equations 



du 
~0t 



— ~L fl (u) = f, 



where L^(u) is described in Section [3J The variational formulation of the problem is given by 

du 



v + a(u,v;fi) = / fv. 

D 01 JD 

In parabolic equations, the basis functions are constructed as in Section [3.21 
For the Galerkin formulation, we write the solution as 

Umsit^lfj) = ^2c t j (t)xi{x)ipj i ' on {x; fx) 

and introduce 

V G n = Sp^( x ^ on )- (37) 
Then the Galerkin formulation is given by 



This framework can also be extended to wave equations. 



G 
on • 



4 Generalization 

The procedure proposed above can be applied to general linear problems. The success of the method will 
depend on the local model reduction that is encoded in a™, s°^, and s™. These bilinear forms need 
to be appropriately defined for a given L^(-). 

One can apply GMsFEM to a linearized nonlinear problem where the operator is frozen at the current 
value of the solution. In this case, one can treat the frozen value of the solution as a scalar parameter on 
a coarse-grid block level. Note that if global model reduction techniques are used, then one will have to 
deal with a very large parameter space and these computations will be prohibitively expensive. 

To demonstrate this concept, we assume that nonlinear equation is linearized 

L,(u n+1 ;u n ) = f. 

For example, for the steady-state Richards' equation, we can consider the following linearization: 

- div(K(x, u n )Vu n+1 ) = f. (38) 

To apply GMsFEM, one can consider u n as a constant parameter, fi — u n , within each coarse-grid 
block. In this case, /i can be regarded as a parameter and we can consider constructing an offline space for 

LM = f. 
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Here fj, £ R is the parameter that represents the average of the solution in each coarse block. In the 
example of the steady-state Richards' equation, u can be assumed to be a constant within a coarse-grid 
block. For general problems (e.g., if the linearization is around Vtx), one needs to use higher dimensional 
parameter space to represent Vu at a coarse-grid level. 

The construction of the offline space will follow the GMsFEM procedure. We will construct a snapshot 
space either by taking local fine-grid functions or by using the value of the parameter corresponding to 
ui that will appear in the linearization of the global system. Furthermore, the offline space is constructed 
via a spectral decomposition of the space of snapshots. Here, appropriate a° ff and s° ff bilinear forms will 
be employed. For example, for Richards' equations, these forms can be chosen as those for linear elliptic 
equation with parameters as described in Section 13.21 In particular, 

<j(^)= E / *(^)|Vx fc | 2 IV>l 2 , 

where a UJi (tp, ip; /x) is the restriction of a{ip,ili;u) onto u>i and Xfe is a partition of unity supported in cok- 
As for s°*, we select 

At the online stage, we consider an approximation of (|38l) with u n replaced by its average in each 
coarse-grid block. We denote this approximate solution by u n+1 

-div(/c(x, (u n ))u n+1 ) = /, 

where (u n ) is the average of u n over the coarse regions. For each parameter value /x, which is the average 
of the solution on a coarse-grid block fi = f u, the multiscale basis functions are computed based on 
the solution of local problem. In particular, for each and for each input parameter and /i = — J u, 
we formulate a quotient to find a subspace of V££ where the space will be constructed for each /x. For the 
construction of the online space, we choose 

a ™('0,V';/Xj) = 51 / n{x,^j)\Vxk\ 2 \^\ 2 , 
k,u lk p l Lj z ^o Juj * 

where Xfc is the partition of unity corresponding to w^, 

= a w*(V>!V>; Mi- 
Note that the value of /x is defined using u n . 

As before, the online space is computed by solving an eigenvalue problem in uji using V^g for the 
current value of u n see (fT5|) . Using dominant eigenvectors, we form a coarse space as in (fT7|) and solve the 
global coupled system following (fT8)l at the current u n . For the numerical example, we will consider the 
coefficient that has an affine representation (see ((H)) which reduces the computational cost associated with 
calculating the stiffness matrix in the online stage (see page Section |3~2")) . 

Remark 8 (Adaptivity in the parameter space). We note that one can use adaptivity in the parameter 
space to avoid computing the offline space for a large range of parameters and compute the offline space 
only for a short range of parameters and update the space. This is, in particular, the case for applications 
where one has a priori knowledge about how the parameter enters into the problem. To demonstrate this 
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concept, we assume that the parameter space A can partitioned into a number of smaller parameter spaces 
Aj, A = IJiAj, where Aj may overlap with each other. Furthermore, the offline spaces are constructed 
for each Aj. In the online stage, depending on the online value of the parameter, we can decide which 
offline space to use. This reduces the computational cost in the online stage. In many applications, e.g., 
in nonlinear problems, one may remain in one of A,; 's for many iterations and thus use the same offline 
space to construct the online space. 

4.1 Nonlinear equation with nonlinearities as a parameter 

In this example, we present an application of GMsFEM to nonlinear equation 

— div(A(a;, u)Vu) = /, 

with u = on dD and 

\{x,u) = \ (x,u) (kiO) +e au K 2 {x)) , (39) 

where k\(x) and Kzix) are defined as in the previous example (see Figure [5]) . The main objective of this 
example is to demonstrate that one can use u in X(x, u) as a scalar parameter within a coarse-grid block. In 
contrast, in global methods, if u in X(x, u) is used as a parameter, it will be a high dimensional parameter. 

We will take the average of u as a parameter in each coarse-grid block as discussed in Section 01 To 
introduce the details of the algorithm, we define the bilinear form a(-, •; •) 

a(u,v;w) — / X(x,w)Vu ■ Vvdx, (40) 
Jd 

and the functional F(-) 

F(v) = / fvdx. (41) 
Jd 

The numerical solution Uh can be approximated to an arbitrary accuracy using a Picard iteration. Starting 
with an initial guess 6 V h , we define the nonlinear fixed point iteration as follows. Given u^J, the next 
approximation is the solution of the linear elliptic equation 

a(u% +1 ,w]u%)=F(w), foredlweV' 1 , (42) 

where u 7 ^ = (u^). This is an approximation of the linear equation — div(fc(x)A(a;, uJJ)Vu^ +1 ) = / with u\ 
being the previous iterate. 

We reformulate the iteration (|4"2l in a matrix form. Define A n by 

a(v, w; uf t ) = w T A n v, for all v, w £ V h , (43) 

and define the vector F by 

F{w) = w T b, for all w G V h . (44) 
Then the equation (|42|) can be rewritten in the following matrix form 

A n ul +1 = b. (45) 

Furthermore, we solve (|45f using GMsFEM with multiscale basis functions computed for all (or a range) 
of in each coarse-grid block. 
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H = 0.1 


N rb = 2 


N rb = 3 


GMsFEM+0 
GMsFEM+1 
GMsFEM+2 


5.079139e- 005(214) 
1.833627e - 005(336) 
1.464074e - 005(460) 


1.759717e- 005(235) 
5.263724e - 006(356) 
4.422483e-006( 477) 



Table 8: Relative errors in energy norm and the coarse space dimension in the last iteration. Here, r] = 10 

In Tabled we present numerical results to study the accuracy of GMsFEM. We take / = 1. We consider 
two cases for N r b, N rb = 2 and N rb — 3. Note that in this case, we have two distinct heterogeneities because 
Ki(x) does not vanish in the linear combination of Ki(x) + e au K2(x). Thus, in this case, N r t = 2 is sufficient 
to capture the heterogeneities across the parameter range. From our numerical results, we observe that 
the errors are small and decrease as we increase the dimension of the local spectral spaces. On the other 
hand, the error is large if MsFEM with one basis function per node is used. This error is 19.21%. 

5 Conclusions 

In this paper, we propose a multiscale framework, the Generalized Multiscale Finite Element Method 
(GMsFEM), for solving PDEs with multiple scales. The main objective is to propose a framework that 
extends MsFEMs to more general problems with complex input space that include parameters, high con- 
trast, and right-hand-sides or boundary conditions. The GMsFEM starts with a family of snapshots for 
the local solutions. These snapshots can usually be generated based on the solutions of local problems 
or simply taking local fine-grid functions. First, based on the space of local snapshots, the offline space 
is constructed. The construction of the offline space involves solving a spectral problem in the space of 
snapshots. This process introduces a prioritization on the space of snapshots across the input space. In 
the online stage of the simulations, for each new parameter and a source term, the online multiscale basis 
functions are constructed efficiently. We discuss various constructions. For example, in the absence of 
parameters, there is no computational work needed in the online stage. When the solution nonlinearly 
depends on the input space parameters, the construction of the online coarse spaces involves solving a 
spectral problem over the offline space. We also discuss the online correction of the reduced solution via 
two-level domain decomposition methods. The optimality of the preconditioners is demonstrated through 
a few examples. We show that the GMsFEM covers some of existing multiscale methods. The general- 
ization of the GMsFEM to nonlinear problems is also considered. We illustrate these methods through 
a few numerical examples. Numerical examples suggest that the proposed framework can be effective in 
studying multiscale problems with large input space dimension and multiple right-hand-sides. 
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